/* This file is part of the Palabos library.
 *
 * The Palabos softare is developed since 2011 by FlowKit-Numeca Group Sarl
 * (Switzerland) and the University of Geneva (Switzerland), which jointly
 * own the IP rights for most of the code base. Since October 2019, the
 * Palabos project is maintained by the University of Geneva and accepts
 * source code contributions from the community.
 *
 * Contact:
 * Jonas Latt
 * Computer Science Department
 * University of Geneva
 * 7 Route de Drize
 * 1227 Carouge, Switzerland
 * jonas.latt@unige.ch
 *
 * The most recent release of Palabos can be downloaded at
 * <https://palabos.unige.ch/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <cstdlib>
#include <iostream>

#include "palabos2D.h"
#include "palabos2D.hh"

using namespace plb;
using namespace std;

typedef double T;

#define DESCRIPTOR descriptors::ShanChenD2Q9Descriptor

template <typename T, template <typename U> class Descriptor>
class RandomInitializer : public BoxProcessingFunctional2D_L<T, Descriptor> {
public:
    RandomInitializer(T rho0_, T maxRho_, Box2D const &boundingBox_, uint32_t const *seed_) :
        rho0(rho0_), maxRho(maxRho_), nY(boundingBox_.getNy()), seed(seed_) {};
    virtual void process(Box2D domain, BlockLattice2D<T, Descriptor> &lattice)
    {
        Dot2D location = lattice.getLocation();
        sitmo::prng_engine eng(*seed);
        plint rng_index = 0;
        for (plint iX = domain.x0; iX <= domain.x1; ++iX) {
            plint globalX = nY * (iX + location.x);
            for (plint iY = domain.y0; iY <= domain.y1; ++iY) {
                plint globalY = iY + location.y + globalX;
                PLB_ASSERT(globalY >= rng_index);
                if (globalY > rng_index) {
                    eng.discard(globalY - rng_index);
                    rng_index = globalY;
                }
                T randNumber = (T)eng() / (T)sitmo::prng_engine::max();
                ++rng_index;

                T rho = rho0 + randNumber * maxRho;
                Array<T, 2> zeroVelocity(0., 0.);

                iniCellAtEquilibrium(lattice.get(iX, iY), rho, zeroVelocity);

                lattice.get(iX, iY).setExternalField(
                    Descriptor<T>::ExternalField::densityBeginsAt,
                    Descriptor<T>::ExternalField::sizeOfDensity, &rho);
                lattice.get(iX, iY).setExternalField(
                    Descriptor<T>::ExternalField::momentumBeginsAt,
                    Descriptor<T>::ExternalField::sizeOfMomentum, &zeroVelocity[0]);
            }
        }
    };
    virtual RandomInitializer<T, Descriptor> *clone() const
    {
        return new RandomInitializer<T, Descriptor>(*this);
    };
    virtual void getTypeOfModification(std::vector<modif::ModifT> &modified) const
    {
        modified[0] = modif::staticVariables;
    };

private:
    T rho0, maxRho;
    plint nY;
    uint32_t const *seed;
};

int main(int argc, char *argv[])
{
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");

    uint32_t seed = 1;

    // For the choice of the parameters G, rho0, and psi0, we refer to the book
    //   Michael C. Sukop and Daniel T. Thorne (2006),
    //   Lattice Boltzmann Modeling; an Introduction for Geoscientists and Engineers.
    //   Springer-Verlag Berlin/Heidelberg.

    const T omega = 1.0;
    const int nx = 400;
    const int ny = 400;
    const T G = -120.0;
#ifndef PLB_REGRESSION
    const int maxIter = 100001;
    const int saveIter = 100;
#else
    const int maxIter = 1001;
#endif
    const int statIter = 100;

    const T rho0 = 200.0;
    const T deltaRho = 1.0;
    const T psi0 = 4.0;

    MultiBlockLattice2D<T, DESCRIPTOR> lattice(
        nx, ny, new ExternalMomentBGKdynamics<T, DESCRIPTOR>(omega));

    lattice.periodicity().toggleAll(true);

    // Use a random initial condition, to activate the phase separation.
    applyProcessingFunctional(
        new RandomInitializer<T, DESCRIPTOR>(rho0, deltaRho, lattice.getBoundingBox(), &seed),
        lattice.getBoundingBox(), lattice);

    // Add the data processor which implements the Shan/Chen interaction potential.
    plint processorLevel = 1;
    integrateProcessingFunctional(
        new ShanChenSingleComponentProcessor2D<T, DESCRIPTOR>(
            G, new interparticlePotential::PsiShanChen94<T>(psi0, rho0)),
        lattice.getBoundingBox(), lattice, processorLevel);

    lattice.initialize();

    pcout << "Starting simulation" << endl;
    for (int iT = 0; iT < maxIter; ++iT) {
        if (iT % statIter == 0) {
            std::unique_ptr<MultiScalarField2D<T> > rho(computeDensity(lattice));
            pcout << iT << ": Average rho fluid one = " << computeAverage(*rho) << endl;
            pcout << "Minimum density: " << computeMin(*rho) << endl;
            pcout << "Maximum density: " << computeMax(*rho) << endl;
        }
#ifndef PLB_REGRESSION
        if (iT % saveIter == 0) {
            ImageWriter<T>("leeloo").writeScaledGif(
                createFileName("rho", iT, 6), *computeDensity(lattice));
        }
#endif

        lattice.collideAndStream();
    }
}
